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High  Resolution  Methods 
for  Time  Dependent  Problems 
with  Piecewise  Smooth  Solutions* 


Eitan  TadrnoE 


Abstract 


A  trademark  of  nonlinear,  time-dependent,  convection-dominated  prob¬ 
lems  is  the  spontaneous  formation  of  non-smooth  macro-scale  features,  like 
shock  discontinuities  and  non-differentiable  kinks,  which  pose  a  challenge  for 
high-resolution  computations.  We  overview  recent  developments  of  modern 
computational  methods  for  the  approximate  solution  of  such  problems.  In 
these  computations,  one  seeks  piecewise  smooth  solutions  which  are  realized 
by  finite  dimensional  projections.  Computational  methods  in  this  context 
can  be  classified  into  two  main  categories,  of  local  and  global  methods.  Local 
methods  are  expressed  in  terms  of  point-values  (-  Hamilton- Jacobi  equations), 
cell  averages  (-  nonlinear  conservation  laws),  or  higher  localized  moments. 
Global  methods  are  expressed  in  terms  of  global  basis  functions. 

High  resolution  central  schemes  will  be  discussed  as  a  prototype  example 
for  local  methods.  The  family  of  central  schemes  offers  high-resolution  “black- 
box-solvers”  to  an  impressive  range  of  such  nonlinear  problems.  The  main  in¬ 
gredients  here  are  detection  of  spurious  extreme  values,  non-oscillatory  recon¬ 
struction  in  the  directions  of  smoothness,  numerical  dissipation  and  quadra¬ 
ture  rules.  Adaptive  spectral  viscosity  will  be  discussed  as  an  example  for 
high-resolution  global  methods.  The  main  ingredients  here  are  detection  of 
edges  in  spectral  data,  separation  of  scales,  adaptive  reconstruction,  and  spec¬ 
tral  viscosity. 

2000  Mathematics  Subject  Classification:  65M06,  65M70,  65M12. 
Keywords  and  Phrases:  Piecewise  smoothness,  high  resolution,  central 
schemes,  edge  detection,  spectral  viscosity. 
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1.  Introduction—  piecewise  smoothness 

A  trademark  of  nonlinear  time-dependent  convection-dominated  problems  is 
the  spontaneous  formation  of  non-smooth  macro-scale  features  which  challenge 
high-resolution  computations.  A  prototype  example  is  the  formation  of  shock  dis¬ 
continuities  in  nonlinear  conservation  laws, 

d 

-u(x,t)  +  Vx-f(u(x,t))  =  0,  U  :=  (iti,  .  .  .  ,Um)T.  (1.1) 

at 

It  is  well  known,  e.g.,  [9],  that  solutions  of  (1.1)  cease  to  be  continuous,  and  (1.1) 
should  be  interpreted  in  a  weak  sense  with  the  derivatives  on  the  left  as  Radon 
measures.  This  requires  clarification.  If  u(x,t)  and  v(x,t)  are  two  admissible  solu¬ 
tions  of  (1.1)  then  the  following  stability  estimate  is  sought  (here  and  below  a,  (3, ... 
stand  for  different  generic  constants), 

||u(-,t)  -u(-,t)||  <  atlK',0)  -u(-,0)||.  (1.2) 

Such  estimates  with  different  norms,  ||  •  ||,  are  playing  a  key  role  in  the  linear 
setting  both  in  theory  and  computations.  For  linear  hyperbolic  systems,  for 
example,  (1.2)  is  responsible  for  the  usual  Instability  theory,  while  the  stability  of 
parabolic  systems  is  often  measured  by  the  L°°-norm,  consult  [14].  But  for  nonlinear 
conservation  laws,  (1.2)  fails  for  any  lA-norm  with  p  >  1.  Indeed,  comparing  u(x,  t ) 
with  any  fixed  translation  of  it,  v(x,  t)  :=  u{x+  h,  t),  the  Lp  version  of  (1.2)  implies 

||A+/tu(-,t)||iP(Rd)  <  at||A+ftu(-,0)||LP(Rd),  A +hu(-,t)  :=  u{ ■  +  h,t)  - 

For  smooth  initial  data,  however,  the  bound  on  the  right  yields  ||  A+hu(-,  t)||ip  < 
cq|/i|,  which  in  turn,  for  p  >  1,  would  lead  to  the  contradiction  that  u{-,t)  must 
remain  continuous.  Therefore,  conservation  laws  cannot  satisfy  the  Lp-stability  es¬ 
timate  (1.2)  after  their  finite  breakdown  time,  except  for  the  case  p  =  1.  The  latter 
leads  to  Bounded  Variation  (BV)  solutions,  ||u(-,t)||w  :=  sup  ||A+/ju(-,  f)||z,1/N  < 
at  <  oo,  whose  derivatives  are  interpreted  as  the  Radon  measures  mentioned  above. 
BV  serves  as  the  standard  regularity  space  for  admissible  solutions  of  (1.1).  A  com¬ 
plete  BV  theory  for  scalar  conservation  laws,  m  =  1,  was  developed  Kruzkov.  Fun¬ 
damental  results  on  BV  solutions  of  one-dimensional  systems,  d  =  1,  were  obtained 
by  P.  Lax,  J.  Glimm,  and  others.  Consult  [2]  for  recent  developments.  Relatively 
little  is  known  for  general  (m  —  1)  x  (d  —  1)  >  0,  but  cf.,  [12]. 

We  argue  that  the  space  of  BV  functions  is  still  too  large  to  describe  the  ap¬ 
proximate  solutions  of  (1.1)  encountered  in  computations.  Indeed,  in  such  compu¬ 
tations  one  does  not  ’faithfully’  realize  arbitrary  BV  functions  but  rather,  piecewise 
smooth  solutions.  We  demonstrate  this  point  in  the  context  of  scalar  approximate 
solutions,  vh(x,t),  depending  on  a  small  computational  scale  h  ~  1/N.  A  typical 
error  estimate  for  such  approximations  reads,  [4] 

II^O^)  —  u(‘)^)llz,i1oc(R)  <  HA-,0)  —  u(-,  0)||Lioc(R)  +ath1^2.  (1.3) 

The  convergence  rate  of  order  1/2  is  a  well  understood  linear  phenomena,  which  is 
observed  in  computations1.  The  situation  in  the  nonlinear  case  is  different.  The 

1  Bernstein  polynomials,  Bjy(u),  provide  a  classical  example  of  first-order  monotone  approxi¬ 
mation  with  Z^-error  of  order  (\\u\\bv /N)1/2 .  The  general  linear  setting  is  similar,  with  improved 
rate  ~  hr^r+1^  for  r- order  schemes. 
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optimal  convergence  rate  for  arbitrary  BV  initial  data  is  still  of  order  one-half, 
[15],  but  actual  computations  exhibit  higher-order  convergence  rate.  The  apparent 
difference  between  theory  and  computations  is  resolved  once  we  take  into  account 
piecewise  smoothness.  We  can  quantify  piecewise  smoothness  in  the  simple  scalar 
convex  case,  where  the  number  of  shock  discontinuities  of  '«(•,  t)  is  bounded  by  the 
finitely  many  inflection  points  of  the  initial  data,  u(x,0).  In  this  case,  the  singular 
support  of  consists  of  finitely  many  points  where  S(t )  =  {x  \  dxu{x,  t )  J.  — oo}. 

Moreover,  the  solution  in  between  those  point  discontinuities  is  as  smooth  as  the 
initial  data  permit  in  the  sense  that  for  arbitrary  I/’s,  L  >>  1,  we  have,  [21] 

sup  \dpu(x,t)\  <  epLT  sup  \dpu(x,0)\+ConstL,  $L(t)  :=  {x  \  dxu(x,t )  >  —  L}. 

x&SL(t)  xGSl(  0) 

If  we  let  d(x,t)  :=  dist(x,S(t))  denote  the  distance  to  S(t),  then  according  to  [19], 
the  following  pointwise  error  estimate  holds,  \vh(x,t )  —  u(x,t)\  <  ath/d(x,t),  and 
integration  yields  the  first-order  convergence  rate 

\\vh(-,t)  -  u(-,t)\\Lioc[R)  <  ath\\og(h)\.  (1.4) 

There  is  no  contradiction  between  the  optimality  of  (1.3)  and  (1.4).  The  former  ap¬ 
plies  to  arbitrary  BV  data,  while  the  latter  is  restricted  to  piecewise  smooth  data  and 
it  is  the  one  encountered  in  actual  computations.  The  general  situation  is  of  course, 
more  complicated,  with  a  host  of  macro-scale  features  which  separate  between  re¬ 
gions  of  smoothness.  Retaining  the  invariant  properties  of  piecewise  smoothness  in 
general  problems  is  a  considerable  challenge  for  high-resolution  methods. 

2.  A  sense  of  direction 

A  computed  approximation  is  a  finite  dimensional  realization  of  an  underlying 
solution  which,  as  we  argue  above,  is  viewed  as  a  piecewise  smooth  solution.  To 
achieve  higher  accuracy,  one  should  extract  more  information  from  the  smooth  parts 
of  the  solution.  Macro-scale  features  of  non-smoothness  like  shock  discontinuities, 
are  identified  here  as  barriers  for  propagation  of  smoothness,  and  stencils  which 
discretize  (1.1)  while  crossing  discontinuities  are  excluded  because  of  spurious  Gibbs’ 
oscillations.  A  high  resolution  scheme  should  sense  the  direction  of  smoothness. 

Another  sense  of  directions  is  dictated  by  the  propagation  of  information  gov¬ 
erned  by  convective  equations.  Discretizations  of  such  equations  fall  into  one  of 
two,  possibly  overlapping  categories.  One  category  of  so-called  upwind  schemes 
consists  of  stencils  which  are  fully  aligned  with  the  local  direction  of  propagating 
waves.  Another  category  of  so-called  central  schemes  consists  of  two-sided  stencils, 
tracing  both  right-going  and  left-going  waves.  A  third  possibility  of  stencils  which 
discretize  (1.1)  ’against  the  wind’  is  excluded  because  of  their  inherent  instability, 
[14].  A  stable  scheme  should  sense  the  direction  of  propagation. 

At  this  stage,  high  resolution  stable  schemes  should  compromise  between  two 
different  sets  of  directions,  where  propagation  and  smoothness  might  disagree.  This 
require  essentially  nonlinear  schemes,  with  stencils  which  adapt  their  sense  of  direc¬ 
tion  according  to  the  computed  data.  We  shall  elaborate  the  details  in  the  context 
of  high-resolution  central  scheme. 
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3.  Central  schemes 

We  start  with  a  quotation  from  [14,  §12.15],  stating  “In  1959,  Godunov  de¬ 
scribed  an  ingenious  method  for  one-dimensional  problems  with  shocks” .  Godunov 
scheme  is  in  the  crossroads  between  the  three  major  types  of  local  discretizations, 
namely,  finite-difference,  finite-volume  and  finite-element  methods.  The  ingenuity 
of  Godunov’s  approach,  in  our  view,  lies  with  the  evolution  of  a  globally  defined 
approximate  solution,  vh(x,tn),  replacing  the  prevailing  approach  at  that  time  of 
an  approximate  solution  which  is  realized  by  its  discrete  gridvalues,  vv(tn).  This 
enables  us  to  pre-process,  to  evolve  and  to  post-process  a  globally  defined  approx¬ 
imation,  vh(x,tn).  The  main  issue  is  how  to  ’manipulate’  such  piecewise  smooth 
approximations  while  preserving  the  desired  non-oscillatory  invariants. 

Godunov  scheme  was  originally  formulated  in  the  context  of  nonlinear  con¬ 
servation  laws,  where  an  approximate  solution  is  realized  in  terms  of  a  first-order 
accurate,  piecewise-constant  approximation 

vh(x,tn)  :=  Ahv{x,tn)  :=  ^vv(tn) \Iv{x),  vv(tn)  :=  f  v(y,tn)dy. 

v  Jl" 

The  cell  averages,  vv,  are  evaluated  over  the  equi-spaced  cells,  /„  :=  {x  \  \x  —  xv  \  < 
h/ 2}  of  uniform  width  h  =  Ax.  More  accurate  Godunov-type  schemes  were  devised 
using  higher-order  piecewise-polynomial  projections.  In  the  case  of  one-dimensional 
equi-spaced  grid,  such  projections  take  the  form 

Vhv(x)  =^2,pv(x)  l/„(a;),  pv(x)  =  vv  +  J  )  + - 

Here,  one  pre-process  the  first-order  cell  averages  in  order  to  reconstruct  accurate 
pointvalues,  vv,  and  say,  couple  of  numerical  derivatives  vj /h,  vj'/h2,  while  the 
original  cell  averages,  {vv},  should  be  preserved,  AhVhVh  =  AhVh.  The  main 
issue  is  extracting  information  in  the  direction  of  smoothness.  For  a  prototype 
example,  let  A+v„  and  A_vv  denote  the  usual  forward  and  backward  differences, 
A±v„  :=  ±(u„±i  —  vv).  Starting  with  the  given  cell  averages,  {u„},  we  set  vv  =  vv, 
and  compute 

vj  =  mm(A+vv,  A_vv),  mm(zi,z2)  :=  sgn^^  s9n(~2)  |^2|}.  (3.1) 

The  resulting  piecewise-linear  approximation  is  a  second-order  accurate,  Total  Vari¬ 
ation  Diminishing  (TVD)  projection,  ||'P/i1’,1(:e)|| w  <  llt’,l(;r)ll w-  This  recipe  of 
so-called  minmod  numerical  derivative,  (3.1),  is  a  representative  for  a  large  library 
of  non-oscillatory,  high-resolution  limiters.  Such  limiters  dictate  discrete  stencils  in 
the  direction  of  smoothness  and  hence,  are  inherently  nonlinear.  Similarly,  nonlinear 
adaptive  stencils  are  used  in  conjunction  with  higher-order  methods.  A  description 
of  the  pioneering  contributions  in  this  direction  by  Boris  &  Book,  A.  Harten,  B. 
van-Leer  and  P.  Roe  can  be  found  in  [10].  The  advantage  of  dealing  with  globally 
defined  approximations  is  the  ability  to  pre-process,  to  post-process  and  in  partic¬ 
ular,  to  evolve  such  approximations.  Let  u(x,t)  =  uh(x,t)  be  the  exact  solution 
of  (1.1)  subject  to  uh(x,tn )  =  VhVh(x,tn).  The  exact  solution  lies  of  course  out¬ 
side  the  finite  computational  space,  but  it  could  be  realized  in  terms  of  its  exact 
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cell  averages,  vh(x,tn+1)  =  ,vv{tn+1)^iu{  x ).  Averaging  is  viewed  here  a  simple 

post-processing.  Two  prototype  examples  are  in  order. 

Integration  of  (1.1)  over  control  volume  /„  x  [tn,tn+1],  forms  a  local  stencil  which 
balances  between  the  new  averages,  {vv(tn+1)},  the  old  ones,  {rv+fc(t”)}>  and  the 
fluxes  across  the  interfaces  along  xv±\/2  x  [tn,tn+1].  In  this  case,  the  solution  along 
these  discontinuous  interfaces  is  resolved  in  terms  of  Riemann  solvers.  Since  one 
employs  here  an  exact  evolution,  the  resulting  Godunov-type  schemes  are  upwind 
scheme.  The  original  Godunov  scheme  based  on  piecewise  constant  projection  is 
the  forerunner  of  all  upwind  schemes.  As  an  alternative  approach,  one  can  realize 
the  solution  u(x,tn+1),  in  terms  of  its  exact  staggered  averages,  {u„+i/2 (fn+1)}. 
Integration  of  (1.1)  over  the  control  volume  Iv+ 1/2  x  [tn,tn+1]  subject  to  piecewise 
quadratic  data  given  at  t  =  tn,  u(x,tn)  =  Vh.vh{x,tn),  yields 

^+i/2(f"+1)  =  -  (yv{tn)  +  vv+i(tnfj  + 

+  l  (vAn  v'v+1(n)  +  ^  (ic++l/2  -  k+1/2)  ,  (3.2) 

where  F"+1^2  stands  for  the  averaged  flux,  F"+1^2  =  J*n  f(u(xv,r)dr/ At.  Thanks 
to  the  staggering  of  the  grids,  one  encounters  smooth  interfaces  xv  x  [tn,tn+1],  and 
the  intricate  (approximate)  Riemann  solvers  are  replaced  by  simpler  quadrature 
rules.  For  second-order  accuracy,  for  example,  we  augment  (3.2)  with  the  mid¬ 
point  quadrature 

K+1/2  =  f(v.(tn+1/2)),  vv{tn+V2)  =  vv{tn)-^-f{vv{tn))'.  (3.3) 

Here,  the  prime  on  the  right  is  understood  in  the  usual  sense  of  numerical  dif¬ 
ferentiation  of  a  gridfunction  -  in  this  case  the  flux  {f(vu(tn))}u.  The  resulting 
second-order  central  scheme  (3. 3), (3. 2)  was  introduced  in  [13].  It  amounts  to  a  sim¬ 
ple  predictor-corrector,  non-oscillatory  high-resolution  Godunov- type  scheme.  For 
systems,  one  implements  numerical  differentiation  for  each  component  separately. 
Discontinuous  edges  are  detected  wherever  cell-averages  form  new  extreme  values, 
so  that  vj(tn)  and  /(^v(fI^)),  vanish,  and  (3. 3), (3. 2)  is  reduced  to  the  forerunner  of 
all  central  schemes  —  the  celebrated  first-order  Lax-Friedrichs  scheme,  [10].  This 
first-order  stencil  is  localized  to  the  neighborhood  of  discontinuities,  and  by  assump¬ 
tion,  there  are  finitely  many  them.  In  between  those  discontinuities,  differentiation 
in  the  direction  of  smoothness  restores  second-order  accuracy.  This  retains  the 
overall  high-resolution  of  the  scheme.  Consult  Figure  1  for  example. 

Similarly,  higher-order  quadrature  rules  can  be  used  in  connection  with  higher- 
order  projections,  [3],  [11].  A  third-order  simulation  is  presented  in  Figure  2.  Finite- 
volume  and  finite-element  extensions  in  several  space  dimensions  are  realized  over 
general,  possibly  unstructured  control  volumes,  x  [tn,  tn+1],  which  are  adapted  to 
handle  general  geometries.  Central  schemes  for  2D  Cartesian  grids  were  introduced 
in  [6],  and  extended  to  unstructured  grids  in  [1],  A  similar  framework  based  on 
triangulated  grids  for  high-resolution  central  approximations  of  Hamilton- Jacobi 
equations  was  described  in  [4]  and  the  references  therein. 

Central  schemes  enjoy  the  advantage  of  simplicity  -  they  are  free  of  (approx¬ 
imate)  Riemann  solvers,  they  do  not  require  dimensional  splitting,  and  they  apply 
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to  arbitrary  flux  functions2  without  specific  references  to  eigen-decompositions,  Ja- 
cobians  etc.  In  this  context,  central  schemes  offer  “black-box  solvers”  for  an  impres¬ 
sive  variety  of  convection-dominated  problems.  At  the  same  time,  central  schemes 
maintain  high-resolution  by  pre  -and  post-processing  in  the  direction  of  smoothness. 
References  to  diverse  applications  such  as  simulations  of  semi-conductors  models, 
relaxation  problems,  geometrical  optics  and  multiphase  computations,  incompress¬ 
ible  flows,  polydisperse  suspensions,  granular  avalanches  MHD  equations  and  more 
can  be  found  in  [16]. 


Figure  1:  Second-order  central  scheme 
simulation  of  semiconductor  device 
governed  by  ID  Euler-Poisson  equa¬ 
tions.  Electron  velocity  in  10'  cm/s 
with  N  =  400  cells. 


Figure  2:  Third-order  central  scheme 
simulation  of  ID  MHD  Riemann  prob¬ 
lem.  with  N  =  400  cells.  The  y- 
magnetic  field  at  t  =  0.2. 


The  numerical  viscosity  present  in  central  schemes  is  of  order  — .  It  is 
suitable  for  the  convective  regime  where  At  ~  Ax,  but  it  is  excessive  when  a  small 
time  step  is  enforced,  e.g.,  due  to  the  presence  of  diffusion  terms.  To  overcome  this 
difficulty,  a  new  family  of  central  schemes  was  introduced  in  [7]  and  was  further 
refined  in  [8].  Here,  the  previous  staggered  control  volumes,  Iv+ 1/2  x  [tn,tn+1], 
is  replaced  by  the  smaller  —  and  hence  less  dissipative,  J„+ 1/2  x  [t",t”+1],  where 
Jv+x/i  '■=  2V+1/2  +  At  x  [a_,a+]  encloses  the  maximal  cone  of  propagation,  a±  = 
au+i/2  =  (min)fe^fe  (/«)■  The  fact  that  the  staggered  grids  are  O(At)  away  from 
each  other,  yields  central  stencils  with  numerical  viscosity  of  order  0(Ax2r~1). 
Being  independent  of  At  enables  us  to  pass  to  the  limit  At  J,  0.  The  resulting  semi¬ 
discrete  high-resolution  central  scheme  reads  vv(t)  =  ~{fv+i/2(t)  —  fv-i/2(t))/ Ax, 
with  a  numerical  flux,  /„_ (-1/2,  expressed  in  terms  of  the  reconstructed  pointvalues, 
V±  =  v±+1/2  =  Vhvh{x„+1/2  ±,t), 


fv+ I/2W  :  — 


a+f(v  )  —  a  f[v+)  , 

— 'A— - iA — L  +  a+a 

a+  —  a 


v+  —  v 
a+  —  a~ 


(3.4) 


2  An  instructive  example  is  provided  by  gasdynamics  equations  with  tabulated  equations  of 
state. 
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Instructive  examples  are  provided  in  Figures  3  and  4. 


Figure  3:  Convection-diffusion  eq.  Figure  4:  Third-order  central  scheme 
ut  +  (u2)x/2  =  (ux/ +  u%.)x  sim-  for  2D  Riemann  problem.  Density  con- 
ulated  by  (3.4), (3.1),  with  400  cells.  tour  lines  with  400  x  400  cells. 

4.  Adaptive  spectral  viscosity  methods 

Godunov-type  methods  are  based  on  zeroth-order  moments  of  (1.1).  In  each 
time  step,  one  evolves  one  piece  of  information  per  spatial  cell  —  the  cell  aver¬ 
age.  Higher  accuracy  is  restored  by  numerical  differentiation  in  the  direction  of 
smoothness.  An  alternative  approach  is  to  compute  higher-order  moments,  where 
the  cell  averages,  vv,  and  say,  couple  of  numerical  derivatives,  vj ,  vj'  are  evolved  in 
time.  Prototype  examples  include  discontinuous  Galerkin  and  streamline-diffusion 
methods,  where  several  local  moments  per  computational  cell  are  evolved  in  time, 
consult  [4]  and  the  references  therein.  As  the  number  of  projected  moments  in¬ 
crease,  so  does  the  size  of  the  computational  stencil.  At  the  limit,  one  arrives  at 
spectral  methods  based  on  global  projections, 

vN(x,t)  =  VNv(x,t)  :=  ^  vk(t)(j>k(x),  %  :=  (v(-,t),<f>k(-))- 

\k\<N 

Here,  <j>k( x)  are  global  basis  functions,  <j>k  —  elkx,Tk(x): ...  and  Vk  are  the  moments 
induced  by  the  appropriately  weighted,  possibly  discrete  inner-product  (•,•).  Such 
global  projections  enjoy  superior  resolution  —  the  error  ||u(-)  —  Vnv(-)\\  decays  as 
fast  as  the  global  smoothness  of  u(-)  permits.  With  piecewise  smooth  solutions, 
however,  we  encounter  first-order  spurious  Gibbs’  oscillations  throughout  the  com¬ 
putational  domain.  As  before,  we  should  be  able  to  pre-  and  post-process  piecewise 
smooth  projections,  recovering  their  high  accuracy  in  the  direction  of  smoothness. 
To  this  end,  local  limiters  are  replaced  by  global  concentration  kernels  of  the  form, 
[5],  KvNvN(x)  =  ^Vk(4>k(x))x,  where  rj(-)  is  an  arbitrary  unit 

mass  G°°[0,1]  function  at  our  disposal.  Detection  of  edges  is  facilitated  by  sepa¬ 
ration  between  the  0(1)  scale  in  the  neighborhoods  of  edges  and  0(hr )  scales  in 
regions  of  smoothness,  KmVn(x)  ~  [u(x)]  +  0(Nd(x))~r .  Here,  [u(a:)]  denotes  the 
amplitude  of  the  jump  discontinuity  at  x  (-  with  vanishing  amplitudes  signaling 
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smoothness),  and  d(x)  is  the  distance  to  the  singular  support  of  v(-).  Two  proto¬ 
type  examples  in  the  27r-periodic  setup  are  in  order.  With  ?;(£)  =  1  one  recovers  a 
first-order  concentration  kernel  due  to  Fejer.  In  [5]  we  introduced  the  concentration 
kernel,  rjexp( £)  exp{/3/£(l  —  £)},  with  exponentially  fast  decay  into  regions  of 

smoothness.  Performing  the  minmod  limitation,  (3.1),  mm{K]^v  n  {x) ,  K^pvn(x)}, 
yields  an  adaptive,  essentially  non-oscillatory  edge  detector  with  enhanced  separa¬ 
tion  of  scales  near  jump  discontinuities,  e.g.,  Figure  7.  Once  macro-scale  features 
of  non-smoothness  were  located,  we  turn  to  reconstruct  the  information  in  the  di¬ 
rection  of  smoothness.  This  could  be  carried  out  either  in  the  physical  space  using 
adaptive  mollifiers,  or  by  nonlinear  adaptive  filters,  op.  In  the  27r-periodic 

case,  for  example,  we  set  'F  *  Vnv(x)  :=  ('S,9,P(- 1  —  -),Vnv(-))  where  'F  is  expressed 
in  terms  of  the  Dirichlet  kernel,  Dp(y)  :=  sl.^fn ^pJy/2)  ’ 

^e,P{y)  =  ^/3(|)-dp((^)’  p  =  pp{y):=ePv2/(y^1)l[-hx\{y)-  (4-1) 

Mollifiers  encountered  in  applications  maintain  their  finite  accuracy  by  localization, 
letting  6  |  0.  Here,  however,  we  seek  superior  accuracy  by  the  process  of  cancellation 
with  increasing  p  j.  To  guarantee  that  the  reconstruction  is  supported  in  the 
direction  of  smoothness,  we  maximize  9  in  terms  of  the  distance  function,  9{x)  = 
d(x),  so  that  we  avoid  crossing  discontinuous  barriers.  Superior  accuracy  is  achieved 
by  the  adaptive  choice  p  d(x)N,  which  yields,  [20], 

\^{d(x),d(x)N}  *  VNv(x)  -  u(x)|  <  Const. {d{x)N)re~')''/d<'x'>N. 

The  remarkable  exponential  recovery  is  due  to  the  Gevrey  regularity  of  pp  £  Q2  and 
Figure  5  demonstrates  such  adaptive  recovery  with  exponential  convergence  rate 
recorded  in  Figure  6.  An  analogous  filtering  procedure  can  be  carried  out  in  the 
dual  Fourier  space,  and  as  before,  it  hinges  on  a  filter  with  an  adaptive  degree  p 

T  *  Tnv(x)  :=  ^2  crp(^jj-^vkelkx,  ap(£)  =  eKP/^2~1\  p  ~  (d(x)N)r,r+1 . 

\k\<N 


Figure  5:  Adaptive  reconstruction  of  Figure  6:  Log  error  for  an  adap- 

the  piecewise  smooth  data  from  its  us-  tive  spectral  mollifier  based  on  N  = 

ing  N  =  128-modes  using  (4.1).  32,  64,  and  128  modes. 
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Equipped  with  this  toolkit  to  process  spectral  projections  of  piecewise  smooth 
data,  we  turn  to  consider  their  time  evolution.  Godunov-type  methods  are  based 
evolution  of  cell  averages.  Cell  averaging  is  dissipative,  but  projections  involving 
higher-order  moments  are  not.  The  following  example,  taken  from  [18]  shows  what 
goes  wrong  with  global  projections.  We  consider  the  Fourier  method  for  27r-periodic 
inviscid  Burgers’  equation,  dtVN{x,  t)+dxSN(v%{-,  t))/ 2  =  0.  Orthogonality  implies 
that  ||t;jv(-,£)||^2  is  conserved  in  time,  and  in  particular,  that  the  Fourier  method 
admits  a  weak  limit,  — *•  v(t).  At  the  same  time,  v(t)  is  not  a  strong  limit, 

for  otherwise  it  will  contradict  the  strict  entropy  dissipation  associated  with  shock 
discontinuities.  Lack  of  strong  convergence  indicates  the  persistence  of  spurious 
dispersive  oscillations,  which  is  due  to  lack  of  entropy  dissipation.  With  this  in  mind, 
we  turn  to  discuss  the  Spectral  Viscosity  (SV)  method,  as  a  framework  to  stabilize 
the  evolution  of  global  projections  without  sacrificing  their  spectral  accuracy.  To 
this  end  one  augments  the  usual  Galerkin  procedure  with  high  frequency  viscosity 
regularization 

^vN(x,t)  +VN^x  •  f(vN{-,t))  =  -N  ^2  (4.2) 

|fc|<JV 

where  ct(-)  is  a  low-pass  filter  satisfying  cr(£)  >  (|£|2s  ^  ■  Observe  that  the 

SV  is  only  activated  on  the  highest  portion  of  the  spectrum,  with  wavenumbers 
|fc|  >  qj\d2s-1)/2s.  Thus,  the  SV  method  can  be  viewed  as  a  compromise  between 
the  stable  viscosity  approximations  corresponding  to  s  =  0  but  restricted  to  first 
order,  and  the  spectrally  accurate  yet  unstable  spectral  method  with  s  =  oo. 


Figure  7:  Enhanced  detection  of  edges  Figure  8:  Legendre  SV  solution  of  in- 
with  v(x )  given  by  v(x)  =  — sgnx  ■  viscid  Burgers’  equation.  Reconstruc- 
cos(x  +  x  ■  sgnx/ 2)l[_7r  jr](a;).  tion  in  the  direction  smoothness. 

The  additional  SV  on  the  right  of  (4.2)  is  small  enough  to  retain  the  formal  spectral 
accuracy  of  the  underlying  spectral  approximation.  At  the  same  time,  SV  is  large 
enough  to  enforce  a  sufficient  amount  of  entropy  dissipation  and  hence,  to  prevent 
the  unstable  spurious  Gibbs’  oscillations.  The  original  approach  was  introduced 
in  [18]  in  conjunction  with  second-degree  dissipation,  s  =  1.  Extensions  to  several 
space  variables,  non-periodic  expansions,  further  developments  of  hyper  SV  methods 
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with  1  <  s  <  oo,  and  various  applications  can  be  found  in  [17].  We  conclude  with 

an  implementation  of  adaptive  SV  method  for  simple  inviscid  Burgers’  in  Figure  8. 
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